A formulation for asphalt concrete air void during service life by adopting a hybrid evolutionary polynomial regression and multi-gene genetic programming

Bitumen, aggregate, and air void (VA) are the three primary ingredients of asphalt concrete. VA changes over time as a function of four factors: traffic loads and repetitions, environmental regimes, compaction, and asphalt mix composition. Due to the high as-constructed VA content of the material, it is expected that VA will reduce over time, causing rutting during initial traffic periods. Eventually, the material will undergo shear flow when it reaches its densest state with optimum aggregate interlock or refusal VA content. Therefore, to ensure the quality of construction, VA in asphalt mixture need to be modeled throughout the service life. This study aims to implement a hybrid evolutionary polynomial regression (EPR) combined with a teaching–learning based optimization (TLBO) algorithm and multi-gene genetic programming (MGGP) to predict the VA percentage of asphalt mixture during the service life. For this purpose, 324 data records of VA were collected from the literature. The variables selected as inputs were original as-constructed VA, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${VA}_{orig}$$\end{document}VAorig (%); mean annual air temperature, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$MAAT$$\end{document}MAAT (°F); original viscosity at 77 °F, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\eta }_{orig,77}$$\end{document}ηorig,77 (Mega-Poises); and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$time$$\end{document}time (months). EPR-TLBO was found to be superior to MGGP and existing empirical models due to the interquartile ranges of absolute error boxes equal to 0.67%. EPR-TLBO had an R2 value of more than 0.90 in both the training and testing phases, and only less than 20% of the records were predicted utilizing this model with more than 20% deviation from the observed values. As determined by the sensitivity analysis, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\eta }_{orig,77}$$\end{document}ηorig,77 is the most significant of the four input variables, while time is the least one. A parametric study showed that regardless of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$MAAT$$\end{document}MAAT, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\eta }_{orig,77},$$\end{document}ηorig,77, of 0.3 Mega-Poises, and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${VA}_{orig}$$\end{document}VAorig above 6% can be ideal for improving the pavement service life. It was also witnessed that with an increase of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$MAAT$$\end{document}MAAT from 37 to 75 °F, the serviceability of asphalt concrete takes 15 months less on average.

Intelligent algorithms, particularly artificial neural networks (ANNs), are increasingly utilized in civil engineering to enhance predictive modeling 28 .Zhao et al. 16 applied multiple machine learning methods, notably finding that support vector regression (SVR) most accurately estimated the initial VA, based on parameters like aggregate gradation and asphalt-aggregate ratio.Similarly, Zavrtanik et al. 1 and Dalhat and Osman 29 leveraged ANNs to analyze asphalt properties and predict specific gravity in asphalt mixes, respectively.Othman 30 also used ANNs to determine key Marshall mix design parameters, highlighting the diverse applications of ANNs in optimizing asphalt mixture designs and improving material property predictions.Due to the black-box nature of ANN and its complex calculations which usually cannot be performed manually, engineers are always in need of new solutions that provide straightforward equations to predict unidentified characteristics and make it easy to predict the outcome of projects by placing new measured values 31,32 .As a result, symbolic regression methods like evolutionary polynomial regression (EPR) and multi-gene genetic programming (MGGP) have become very prominent and widely employed in recent years.
An extended regression method called EPR incorporates the advantages of both least square regression and genetic programming (GP) regression approaches 33 .EPR has effectively modeled various engineering challenges, such as the behavior of unsaturated soils, slope stability, soil permeability, asphalt flow number, and mechanical properties of stabilized sandy soil [34][35][36][37][38][39][40][41] .In this study, the teaching-learning based optimization (TLBO) algorithm, a population-based evolutionary method that emulates the teacher-student interaction, was used to refine the constants in EPR models, enhancing their accuracy and efficiency without the need for tuning parameters other than iteration count and population size 42,43 .Additionally, MGGP has been utilized as a machine learning technique to generate nonlinear predictive formulas efficiently.MGGP, founded on principles of genetics and natural selection by Koza [45][46][47] , leverages parse trees instead of fixed-length binary strings, providing the flexibility needed for complex problem-solving across various domains without being bound by prior data patterns 44,48 .This approach has addressed numerous engineering challenges effectively [49][50][51][52][53][54][55] .
By considering what has been briefly explained above, there is a gap in the development of models that can provide a simple and practical equation for predicting aged VA values more accurately, which can facilitate the www.nature.com/scientificreports/design process of asphalt concrete, enhance its performance prediction, and finally, improve the safety of the roads.It is the goal of this study to address these gaps by: (I) compiling a comprehensive database of asphalt concrete VA during the service life utilizing published sources, (II) implementing hybridized EPR by applying the TLBO algorithm and MGGP for modeling VA, (III) comparing the models presented in this study to previously developed models, and (IV) identifying the importance of input variables on VA of asphalt mixture during its service life and conducting parametric study based on each inputs variation in a specific range and also predicting VA.Through this research, a further understanding of VA during the service life will be gained and it will pave the way for a more accurate simulation of asphalt concrete performance during the service life.

Methodological background EPR
Over the past two decades, advances in computational tools have considerably facilitated the use of computational intelligence techniques in fields like data mining and soft computing 33 .Among these techniques, EPR stands out due to its integration of traditional regression and GP elements 56 , making it effective for modeling complex environmental phenomena 38,[57][58][59][60] .First introduced by Giustolisi and Savic 61 , EPR operates through a two-step process that constructs symbolic models by determining structure and refining relationships between inputs and outputs 62 .This method evaluates the accuracy of its output functions through the correlation coefficient and adjusts its models based on the desired accuracy and other criteria like the number of generations or sentences in the model 61 .The EPR method's efficacy lies in its ability to produce high-fit models quickly, employing polynomial relations with a limited number of sentences to avoid unnecessary complexity while maintaining accuracy 33 .Various forms of the equation are used to fit the data, involving different powers and user-defined transformations of the input variables to best model the relationships within the data.Each version of the equation demonstrates the flexibility of EPR in accommodating different data structures and analysis needs, with parameters adjusted to refine predictions and improve model fit 63 .This adaptive, iterative method allows EPR to remain robust across various applications, dynamically adjusting to new data and findings.Figure 1 displays the flowchart of the EPR method.

Start
Initialize input matrix

Create initial population of exponent vectors randomly
Assign exponent vectors to corresponding columns of input matrix to create a population of mathematical structures Evaluate coefficients using least square method to create a population of equations

Output results End
Evaluate fitness of equations in the population

Select individuals from mating pool of exponent vectors
Select two exponent vectors to conduct crossover Select one exponent vector to perform mutation

GA Tool
Is termination criterion satisfied?

TLBO
The TLBO algorithm, developed by Rao et al. 65 , is inspired by the educational influence of teachers on students in a classroom setting.TLBO operates in two main phases as the teacher phase and the student phase.In the teacher phase, the algorithm simulates a classroom where one "teacher" (the best solution) aims to improve the overall learning (performance) of a group of "students" (solutions).Each student's performance on various "subjects" (design variables) is evaluated, and the best performing student's results influence the rest of the population.The key mathematical representation of this phase involves adjusting each student's solution based on the difference between the best solution and the class average, modulated by a random factor (r i ) and a teaching factor (T F ), which determines the magnitude of learning: where, X j,kbest,i is the result of the best student in subject j.The teaching factor T F specifies how the average should be changed, and r i contains a random number ranging from 0 to 1. Also, X'j,k,i is the modified value of X j,k,i .When X j,k,i ′ provides a superior function value, it is acceptable.The student phase follows, where each student learns from peers, further refining solutions through the mutual interaction.This interaction encourages diversity in the solution pool, preventing premature convergence to local optima.Students randomly exchange information, adjusting their solutions based on those who have performed better: In which X ′′ j,P,i is valid as long as resulting in a superior function value 66 .TLBO is noted for not requiring specific algorithmic parameters like crossover rates or mutation probabilities, typically necessary in other evolutionary algorithms.This simplification, along with its dual-phase learning approach, makes TLBO a unique and effective tool for solving optimization problems across various engineering disciplines 66 .The algorithm's effectiveness lies in its ability to simulate the dynamic and interactive learning environment of a classroom, continuously enhancing the solution quality through collaborative and competitive processes.The TLBO algorithm is illustrated in Fig. 2.

MGGP
GP, devised by Koza 45 , employs the principles of genetics and natural selection to effectively solve problems across various domains 46,47 .This method uses parse trees, flexible structures that replace fixed-length binary strings, allowing it to adapt to different problem types without being restricted by domain-specific patterns 45,48 .GP operates through a series of evolutionary processes including selection, crossover, and mutation, to optimize solutions based on fitness, which is measured against a predefined error metric 48,67,68 .The flexibility of GP is enhanced by its ability to handle a variety of input variables and functions, termed terminals and primitive functions, respectively, which guide the development of solutions 53 .Moreover, the GP's effectiveness in producing reliable and diverse solutions is achieved by managing its evolutionary operators, with a particular emphasis on maintaining genetic diversity through mutation to prevent premature convergence 48 .MGGP extends GP by integrating multiple parse trees into a cohesive model, where each tree represents a distinct "gene" contributing to the model's output 69 .This approach not only allows for the modeling of complex, non-linear relationships but also ensures that the model remains interpretable by maintaining a balance between complexity and parsimony.The MGGP's ability to combine linear regression with non-linear transformations helps capture intricate interactions within data, providing a robust framework for predictive modeling that is both accurate and efficient 70,71 .This blend of GP and advanced regression techniques makes GP and MGGP powerful tools in the field of computational intelligence, capable of tackling a wide range of engineering and data science problems.Figure 3 shows a typical flowchart for an MGGP procedure.

Data collection
The database utilized in this research is from a report by the Transportation Research Board (TRB) of the National Research Council 26 .Among the three parts of the data used in this study, the first comes from GPS-LTPP database, which includes 72 data, the second, contained 177 data collected by Mirza 27 , and the third, derived from NCAT database 26 , which includes 75 data.Further information regarding each of the data sets can be found in the mentioned literature, including project location details, binder grades, or number of sections on each test road.Additionally, in this study, the aged VA (%) was selected as the output of the model, and original as-constructed VA, VA orig (%); mean annual air temperature, MAAT (°F); original viscosity at 77 °F, η orig,77 (Mega-Poises); and time (months) were considered as inputs.VA orig provides a crucial baseline, allowing for the quantification of changes in VA structure due to traffic and environmental conditions.MAAT is vital for predicting thermal effects on asphalt properties, as temperature significantly impacts the viscoelastic behavior of asphalt.η orig,77 allows for the assessment of the asphalt's flow characteristics at standard temperatures, influencing how it responds to stress over time.Tracking time enables modeling the evolution of these properties, leading to precise predictions (4) Difference_Mean j,k,i = r i X j,kbest,i − T F M j,i (5) X j,k,i ′ = X j,k,i + Difference_Mean j,k,i (6) X ′′ j,P,i = X j,P,i ′ + r i X j,P,i ′ − X j,Q,i ′ , ifX total,P,i ′ < X total,Q,i ′ (7) X ′′ j,P,i = X j,P,i ′ + r i X j,Q,i ′ − X j,P,i ′ , ifX total,Q,i ′ < X total,P,i ′ Vol.:(0123456789)  4 also depicts the histogram and cumulative frequency for input and output parameters.It can be seen that the amount of VA orig is between 3.70 and 20.69%, time varies between 0 and 288 months, η orig,77 falls between 0.12 and 11.97 Mega-Poises, and MAAT ranges from 37.00 to 74.90 °F.Also, the measured VA during the service life changes in the range of 1.14 and 20.40%.A variety of changes in each of the inputs and outputs indicates that the database is diverse, and the models developed over it are generalizable 72 .It should be noted that the original VA depend on the mix design and compaction during construction.Higher initial VA leads to a more noticeable fluctuation in VA compared to lower initial VA.Over time, under various conditions like mix type, traffic, and environment, a final VA value is reached 27 .This final VA is often referred to as "voids at refusal".Heitzman et al. 73 demonstrated that the asconstructed or original VA in selected LTPP sections alters between 1% to more than 14%.

Model performance evaluation criteria
Selecting appropriate error measures like R 2 , RMSE, MAE, MAPE, OBJ, and a20-index for evaluating the accuracy and reliability of models in this study was driven by their specific relevance to the objectives and data characteristics 74,75 .R 2 is used to determine the proportion of variance in the dependent variable that is predictable from the independent variables, making it essential for assessing the model fit.RMSE and MAE provide insights into the average model prediction errors, with RMSE being more sensitive to outliers, thus offering a robust Initialize number of students (population), termination criterion

Calculate mean of each design variables
Identify best solution (teacher)

Modify solution based on best solution X new = X old + r . (X teacher -(T F )mean)
Is new solution better than existing?Accept Reject Select any two solutions randomly X i and X j Is X i better than X j ?
Is new solution better than existing?Accept Reject Is termination criteria satisfied?measure of error magnitude 76 .MAPE aids in understanding error in terms of percentage, which is particularly useful for comparisons across different scales or units.The OBJ function combines several statistical measures to provide a comprehensive error summary 77 , while the a20-index offers a unique perspective by focusing on   www.nature.com/scientificreports/ the accuracy within the 20th percentile 78 , crucial for scenarios where smaller prediction errors are critical.Each of these metrics was chosen to ensure a holistic evaluation of the model's performance, directly supporting the study's aims to develop reliable predictive models for real-world applications.These statistical indicators are summarized below: in which, Y obs and Y pre indicate the observed and predicted values, respectively, where the bar signs over the above-mentioned variables represent the average value; N is the number of records; m20 is the number of records in which their Y obs /Y pre ratio falls between 0.80 and 1.20; and tr and tst refer to the training and testing datasets, respectively.

EPR modeling
The EPR method was used to model the data after it was split into training (75%) and testing (25%) data.This process involved internal functions, including logarithmic, exponential, hyperbolic tangent, and hyperbolic secant.Based on the minimum modeling error, the best model was selected for each function and structure.Also, the presentation type of EPR models was considered in hyperbolic secant, hyperbolic tangent, and logarithmic models as Eq. ( 14), and in exponential models as Eq. ( 15).Table 2 outlines other details set for the development of each model.
The representative range mentioned in Table 2 determines the linearity or non-linearity of the model as well as the direct or inverse dependence between input parameters and output parameters.
According to Eqs. (16-19), the optimal models for different functions can be found as follows:  In accordance with the training and testing data, Fig. 5a-d demonstrate each EPR model's ability to predict VA based on measured and predicted values.In a chart, a higher R 2 value reveals a more efficient model and fewer scattering points than a line of 100% agreement.In comparison with other models, the hyperbolic secant model is selected as the best model based on its highest R 2 and lowest RMSE values.In the case of the hyperbolic secant model, which was chosen as the best model, the R 2 values for the training and testing data are 0.8956 and 0.8995, respectively.Also, the RMSE values of 1.0042 and 1.0288, respectively, for the training and testing data, point out an acceptable level of accuracy for the model.An a20-index has been introduced as a new parameter of physical engineering that refers to the number of samples where the difference between the observed and predicted values is less than 20% 79 .Based on Fig. 5, the a20-index equals 0.844 and 0.815 for the hyperbolic secant model in the training and testing phases, respectively.Moreover, in both the training and testing phases, the exponential model gave the worst performance based on a20-index values of 0.65 and 0.61, respectively.

Optimizing EPR model with TLBO algorithm
In this study, the TLBO algorithm was employed to improve the validity of the hyperbolic-secant-derived EPR model by fine-tuning constant coefficients, assuming the form of the model is fixed.Accordingly, the objective ( 16)  A number of statistical parameters were evaluated for the overall dataset to assess the precision of the EPR-TLBO model in making forecasts.The results for different statistical parameters are listed in Table 3. Comparing the R 2 values of EPR-TLBO and hyperbolic secant models reveals that the first model is 1.9% more superior than the second in the training phase, but this superiority drops to 0.78% in the testing phase.Furthermore, RMSEs for the EPR-TLBO and hyperbolic secant models are 0.989 and 1.029%, respectively.MAEs of the EPR-TLBO and hyperbolic secant models are 0.685% and 0.602% in the training phase, respectively, whereas their MAPEs are, respectively, 9.566% and 9.593% in the testing phase.

MGGP modeling
Genetic operators, number of populations, tournament size, number of genes, ERC, crossover, mutation probabilities, and elite fraction are the primary MGGP factors whose rates vary and play an important role in the MGGP algorithm 53,80 .Researchers have proposed 31,81 that MGGP hyper-parameters must be tuned in a trialand-error manner to build a VA equation based on VA orig , MAAT , η orig,77 , and time .An optimum setting of MGGP hyper-parameters is obtained through a trial-and-error procedure and by taking into account previous research findings.
As inputs to the model, terminal sets are unrelated variables.As a result, we used trigonometric and arithmetic functions near our problem in the present study.Moreover, the functions of + , − , × , /, sqrt, exp, Ln, ^2, ^3, 3Rt, sin, cos, and tanh were utilized to design a superior MGGP model.It is necessary to link genes in the multi-genic chromosomes, as well as create large and elaborate gene trees.Here, addition operators (+) are (time + 0.0001) www.nature.com/scientificreports/employed to link genes, as these operators give better results than the rest (e.g.− , × , /).Trial-and-error can be used to obtain other MGGP parameters, such as the number of populations, generations, and runs.In the study, 100 populations, 250 generations, and five runs were considered.By changing the number of genes and keeping the other settings parameters fixed, the trial-and-error procedure was followed to select the optimal architecture for the MGGP model.Additionally, the number of nodes was fixed at "Inf" so that the MGGP algorithm had more freedom to build models.Table 4 outlines the performance of various MGGP models with different numbers of genes from 1 to 19.As can be seen, utilizing 15 genes yields the highest testing R 2 value (0.897).Additionally, the 15-gene model exhibited lower errors in terms of MAE, RMSE, and MAPE values compared to the other models.As shown in Fig. 6, the best chromosome is a tree structure made up of 15 genes connected by addition functions (+).Finally, VA can be predicted with Eq. ( 21).Also, the predicted VA values of the MGGP model and measured VA of the training and testing datasets can be seen in Fig. 7. Hence, MGGP is fairly accurate with a20-index values of 0.84 and 0.815 for the training and testing phases, respectively.

VA orig Sqrt
By combining different statistical indicators for the training and testing datasets, OBJ allows for assessing the performance of a model 82 .When an OBJ high value occurs, it means that a model performs poorly compared to other models 83 .According to Fig. 8, the EPR-TLBO model performs best with an OBJ of 0.81.With an OBJ value of 0.889, the hyperbolic secant model is the next priority, followed by the MGGP model with 0.891.As a further measure of the validity of the models, the scatter index (SI) as well as the Nash-Sutcliffe efficiency (NSE) coefficient were used: Observed values are represented by Y obs and predicted values are denoted by Y pre , while the bar symbols over the respective values illustrate the average value, and N is the total number of records; When the SI value of a model is less than 0.1 or the NSE value is greater than 0.7, the prediction capability of the model is excellent.A model can be considered good if it can achieve SI values between 0.1 and 0.2 or NSE values between 0.65 and 0.75, while if the SI value is between 0.2 and 0.3, or if the NSE value is between 0.5 and 0.65, it can be considered fair, and a model with lower SI and NSE values will perform poorly 72,74 .Based on Fig. 9, the NCHRP, EPR-TLBO, and MGGP models have NSE values above 0.75, indicating excellent learning and generalization accuracy for  www.nature.com/scientificreports/ the NCHRP, EPR-TLBO, and MGGP models.Meanwhile, the NCHRP, EPR-TLBO, and MGGP models have SI values between 0.1 and 0.2, which demonstrates that they are good predictors.The results also reveal that the EPR-TLBO model is more accurate due to having a lower value of SI and a higher value of NSE.

M ir z a N C A T N C H R P E P R -T L B O M G G P
The performance of each model was compared based on the Taylor diagram displayed in Fig. 10.RMSE, R 2 , and STD are statistical parameters used to compare predictions with observations.In the plot, the standard deviation is represented as a circular line connecting the vertical and horizontal axes.Moreover, the horizontal circular green dots demonstrate the RMSE values while the radial blue lines exhibit the R 2 values.A comparison of the EPR-TLBO model to the observed data reveals that it performs the best.It matches the observed data the closest among all the models.
In accordance with Fig. 11, the logarithms of the predicted to measured records ratios for all the VA models for the testing data and their normal density functions were analyzed.Density functions with a higher amplitude show higher accuracy 82 .As can be seen, the maximum density is found close to the zero value of the x-axis for the EPR-TLBO and MGGP models, respectively, highlighting their high testing efficiency.

Sensitivity analysis
Validating and testing the robustness of the model for unknown data require different analyses of the proposed model.By examining the sensitivity of the input variables, we can determine how the output of a model is affected by the input variables.Here, the sensitivity analysis was conducted according to the method proposed by Gandomi et al. 84 and Javed et al. 85 .In this method, input variables are calculated as a percentage of the predicted VA www.nature.com/scientificreports/through the sensitivity analysis.This calculation is carried out with Eqs.(24 and 25), which are used to determine the sensitivity percentage of VA for any input variable: where f min p i and f max p i are the lowest and highest values of the outputs, respectively, calculated from the i th input variable, assuming that all other input variables are held constant in their mean values.In SP = 0, the amount of contribution is the lowest, while SP = 100 represents the highest contribution.Figure 13 features the results of the sensitivity analysis of the input variables for the EPR-TLBO model.As can be seen, VA orig (%) and MAAT (ºF) provide the largest and smallest contributions to the asphalt concrete mixture VA, respectively.Overall, the importance of features can be prioritized as VA orig (%) > time (months) > η orig,77 (Mega-Poises)> MAAT (ºF).This finding is consistent with the results of the sensitivity analysis for the NCHRP model 26 .The graphical sensitivity analysis employing the NCHRP model revealed that two parameters VA orig and time strongly affect the VA value, and the MAAT parameter has a much smaller effect on VA during the service life 26 .

Parametric study
The life expectancy of conventional road pavement is between 20 and 30 years 86 .Asphalt concrete pavements are commonly damaged by fatigue, along with thermal cracking and rutting 87 .The longitudinal depression along the wheel path, called rutting, usually has upheavals on both sides 88 .Due to the high as-constructed VA content of the material, it is expected that VA will reduce over time, causing rutting during initial traffic periods 88 .Eventually, the material will undergo shear flow when it reaches its densest state with optimum aggregate interlock or refusal VA content, typically between 2 and 3% 89 .Moreover, there should be enough flexibility in the asphalt to resist distress.When asphalt-coated aggregate particles are compacted, they become stable and resistant to various types of deformation while improving the mixture's durability and reducing its permeability 90 .Several researchers concluded that dense-graded asphalt mixes should not contain less than 3.0% VA during the service (24) life due to subsequent pavement distresses [91][92][93][94] .Based on the previous studies, adding just 1.0% more VA to the asphalt mix may result in a 35% decreased pavement fatigue life and doubled permeability 95 .Also, an increase of 1.0% in VA over the base level of 7.0% tends to reduce the pavement life by 10% 96 .Thus, engineers and road designers can make accurate and cost-effective decisions by having a fairly accurate estimate of the service life of asphalt pavement in exchange for different influencing variables such as VA orig (%), η orig,77 (Mega-Poises), and MAAT (°F).However, researchers are reluctant to conduct experimental studies because of time constraints and the cost of preparing numerous asphalt core samples and doing lab tests.With the development of prediction models, parametric analysis can also be performed along with determining how input variables affect output.To achieve this, two parameters of time and MAAT were considered as being between the minimum and maximum levels in Fig. 14, and other parameters such as η orig,77 of 0.1, 1.0, and 3.0 Mega-Poises, and VA orig of 6, 9, and 12% were assumed as constants.After that, VA was calculated according to the desired changes in parameters.
In a general view, the VA percentage is negatively correlated to time and MAAT , whereas VA orig and η orig,77 are directly related to the VA amounts.To find an expectation for the service life of asphalt concrete, a plane with a constant VA value of 3.0 cut the main mesh surface.As a result, the intersection of these surfaces determines the approximate service life of asphalt concrete.In Table 6, low and high MAAT s of 37 and 75 °F are compared with their respective calculated service life.At η orig,77 of 0.1 Mega-Poises, increasing VA orig from 6 to 9 and 12% adds 70 and 140 months to the asphalt concrete service life equally at the minimum and maximum MAAT s, respectively.Also, at η orig,77 = 1.0 Mega-Poises, a 1.5-fold increase in VA orig extends the service life by 2.8 to 3.1 times at MAAT of 37 and 75 °F, respectively.Meanwhile, when VA orig doubles, at different MAAT s, the pavement serviceability takes 4.6-5.4times longer.In the case of η orig,77 = 3.0 Mega-Poises, for VA orig higher than 6%, the study was unable to provide a pavement life expectancy estimate.Table 6 enables a nuanced analysis, allowing designers to assess how various factors might influence the longevity and durability of asphalt under typical service conditions.By integrating these parameters into their evaluations, designers can make more informed decisions regarding material choices and construction strategies, ultimately enhancing the performance and sustainability of pavement installations.In the analysis of the service life, only the rutting life and bleeding of asphalt mixture have been considered, and other failures such as fatigue damage, moisture sensitivity, and aging due to high VA have not been taken into account.

Limitations
The current limitations in the data and methodology for the pavement performance modeling are largely due to insufficient comprehensive datasets, particularly lacking in detailed traffic load and environmental data.To enhance future research, there is a need for more robust data collection systems, incorporating real-time traffic monitoring and environmental conditions, alongside more frequent pavement condition assessments.Additionally, integrating advanced machine learning techniques like deep learning, adding more predictive variables such as material properties and maintenance history, and employing rigorous model validation strategies can remarkably improve the accuracy of these models.By expanding the depth and breadth of data and employing crossdisciplinary approaches, future studies can develop more precise and reliable predictive models, thus enhancing decision-making in pavement management and extending the service life of pavement infrastructures.In order to specify the limitations of the dataset used in this study, the database specifications including the number of data points and the range of independent variables are listed in Table 7.

Conclusions
This article aimed to predict VA values of asphalt concrete mixtures throughout their service life by adopting a hybrid EPR approach with the TLBO algorithm and GP.As part of this effort, 324 data from the literature have been collected regarding VA during the service life.As input variables, we selected VA orig (%), time (months), η orig,77 (Mega-Poises), and MAAT (°F).Below is a summary of the findings: • TLBO has shown to be capable of optimizing the EPR coefficients in an efficient manner.
• In the testing and training phases, the EPR-TLBO model is superior to the GP model according to statistical indicators.All of the proposed models in this study were correctly trained and all of their predicted values were highly correlated with low error rates, based on a comparison between the EPR-TLBO and GP models with those from the literature that was conducted in this study.• The EPR-TLBO and GP models were able to predict over 80% of records with a maximum deviation of 20% from the actual results.• The newly developed models exhibited excellent performance with SI between 0.1 and 0. • The EPR-TLBO models have the superior capability when compared to other existing models because of the interquartile ranges of absolute error boxes equal to 0.67%.• The highest normal density function values were found near the zero point of the logarithms of the predicted to measured results ratios for both the EPR-TLBO and GP models, respectively, highlighting their remarkable testing efficiency.
• A parametric study revealed that regardless of MAAT , η orig,77 of 0.3 Mega-Poises and VA orig above 6.0%can be ideal for improving the pavement service life.It was also observed that with an increase of MAAT from 37 to 75 °F, the serviceability of asphalt concrete takes 15 months less on average.• EPR-TLBO, due to its robust optimization capabilities, may excel in complex scenarios but can be prone to overfitting if not properly regulated.MGGP, offers a flexible model structure that can adaptively fit diverse data patterns, but it may also become overly complex, making it difficult to interpret.Traditional regression models, while generally easier to understand and implement, might lack the nuanced fitting capabilities of more advanced algorithms, potentially leading to under fitting in complex datasets.Each model's performance thus reflects a trade-off between accuracy and model interpretability, highlighting the importance of selecting the right model based on specific project requirements and data characteristics.• The new models proposed in this study allow for more accurate predicting of VA of asphalt concrete during the service life.More accurate predicting of VA during the service life results in more accurate predicting of dynamic modulus change during the service life which is crucial for the pavement analysis and design for both reconstruction and rehabilitation projects.Also, accurate models for predicting VA of asphalt concrete during the service life are beneficial to detect the premature damage of asphalt concrete layers related to low VA such as bleeding and rutting.

Figure 3 .
Figure 3.A typical flowchart for an MGGP procedure.

Figure 9 .
Figure 9.Comparison of SI and NSE calculated values.

Figure 12 .
Figure 12.Absolute error boxes of various EL methods.

Figure 14 .
Figure 14.Analysis of impact of input variables on aged VA of asphalt concrete mixtures.
2 and NSE greater than 0.75.• A comparison of the estimated values and their correlation with the observed values suggests that the predic- tions of the EPR-TLBO model are closer to those of Mirza, NCAT, and NCHRP.A high OBJ value indicates poor existing model's performance over those proposed in this article.• As determined by the sensitivity analysis, η orig,77 (Mega-Poises) is the most significant of the four input vari- ables.

Table 1 .
Descriptive statistics of training and testing data.

Table 2 .
Details of tuning parameters for models.

Table 3 .
Accuracy and performance of each EPR model.

Table 4 .
A comparison of various MGGP models for training and testing records.Significant values are in bold.

Table 5 .
Statistical values of developed models.Significant values are in bold.

Table 5
compares the performance of the models proposed in this article with the equations available in the research literature for predicting VA values across all 324 data records.Generally, both developed models in this study were superior to those developed by Mirza, NCAT, and NCHRP.As demonstrated, the best model to predict VA of asphalt concrete is EPR-TLBO model with the highest R 2 (0.91) and lowest RMSE (0.936%), followed by the MGGP model with an R 2 value of 0.897 and RMSE value of 1.0%; NCHRP with an R 2 value of 0.865 and an RMSE value of 1.114%; Mirza with an R 2 value of 0.8098 and an RMSE value of 1.413%; and NCAT with an R 2 value of 0.789 and an RMSE value of 1.454%.Calculated MAE and MAPE values signify that the accuracy of the models is prioritized as follows: EPR-TLBO > MGGP > NCHRP > Mirza > NCAT.It should be noted that the performance of the MGGP model in predicting VA values is not far from the EPR-TLBO model.

Table 6 .
An overview of parametric study′s results.

Table 7 .
Characteristics of dataset for developing models in this study.